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Abstract 

We study the propagation of small amplitude waves superimposed on a 
large static deformation in a nonlinear viscoelastic material of differential 
type. We use bulk waves and surface waves to address the questions of 
dissipation and of material and geometric stability. In particular, the 
analysis provides bounds on the constitutive parameters and on the pre- 
deformation that ensure linearized stability in the neighbourhood of a 
large pre-stretch. This type of result is relevant to the imaging of biological 
soft tissues using acoustical techniques, where pre-deformation is known 
to increase contrast and reduce de-correlation noise. 

1 Introduction 

In many important technological applications, polymeric materials — such as 
the elastomers used in engine mounts or bridge bearings — are subject to large 
deformations, and infinitesimal theories are not suitable for modelling their me- 
chanical response. This is true also for complex biomaterials 'in service' such 
as ligaments, tendons, skin, arteries, and other biological soft tissues that have 
several mechanical features in common with elastomeric polymers. 

An adequate modelling of rubber-like materials and of biological soft tis- 
sues subject to large deformations requires the use of the theories of nonlinear 
elasticity and nonlinear viscoelasticity. Whereas the theory of nonlinear elas- 
ticity is a well-developed chapter of solid mechanics, the theory of nonlinear 
viscoelasticity is still in its infancy. Relatively few studies have been carried out 
beyond establishing basic constitutive characterizations and their general ther- 
modynamical implications. In particular, there is a paucity of complete studies 
of the propagation of mechanical waves. Beyond the literature dedicated to ac- 
celeration waves and universal motions, we find few papers dedicated to finite 
amplitude waves and to small amplitude waves superimposed on finite defor- 
mations in viscoelastic solids (of course, the situation is different for waves in 
viscoelastic fluids). 

Antman and Seidman pQ provided a detailed mathematical study of large 
shearing motions of nonlinear ly viscoelastic slabs (see also Rajagopal and Sac- 
comandi [27j for some exact solutions in a similar framework). Recently, Hayes 



and Saccomandi [HJ [331 H3] and Destrade and Saccomandi [HI M HD] obtained 
some results for finite amplitude motions and waves in some classes of nonlin- 
ear viscoelastic materials. Earlier, Hayes and Rivlin [T71 [TH1 [HO HO] had 
established some general results for the theory of small motions superposed on 
a large deformation in nonlinear viscoelastic solids (see also a recent note by 
Saccomandi [5H] concerning such waves in a special class of materials). 

The situation is, of course, completely different in the linear theory of vis- 
coelasticity where the study of bulk and surface waves is a well-developed subject 
with a wealth of results obtained over the years. However, this linear framework 
does not meet the needs of the actual technological advances in non-invasive 
techniques of investigation and medical imaging. These techniques, based on 
ultrasound [30j , are bringing wave motion to the forefront of imaging and ther- 
apy in many areas of medicine. At the same time, the apparatus used must 
now rely on nonlinear constitutive assumptions in order to account for the pre- 
loading and large stretches found in living soft tissues; see the recent review by 
Hoskins [24] . In particular, we emphasize that pre-loads and pre- deformations 
are fundamental for reducing the dynamic range of object stiffness. Indeed, 
compression of soft tissues before imaging increases contrast and reduces de- 
correlation noise [HI [13] . 

Here we study the propagation of small amplitude waves in certain isotropic 
and incompressible nonlinear ly viscoelastic solids, with a view to investigating 
their stability when subject to large deformations. The solids under consid- 
eration are characterized by a Cauchy stress tensor T depending only on the 
Cauchy-Green deformation tensor B and on the symmetric part D of the ve- 
locity gradient. This class of materials is usually referred to as simple materials 
of differential type of grade 1 ; see Truesdell and Noll [31] . This is a basic class 
of models in nonlinear viscoelasticity; it accounts for classical effects like creep 
and recovery, as in Kelvin- Voigt linear viscoelasticity, but cannot describe stress 
relaxation. This class contains the so-called Mooney-Rivlin viscoelastic mate- 
rial [5J and the incompressible version of the model proposed by Landau and 
Lifschitz in their book on the theory of elasticity [25] . 

In Section 2 we summarize the basic governing equations and constitutive 
assumptions for these materials. We devote Section 3 to deriving the general 
form of the incremental equations of motion in a deformed viscoelastic solid. 
Then we specialize the analysis to two-dimensional motions and find conditions 
for time-averaged dissipation in time-periodic homogeneous motions. In Sec- 
tion 4 we study bulk wave propagation and material stability; we find that the 
combination of time-averaged dissipation and strong ellipticity of the static de- 
formation results in material stability. In Section 5 we consider surfaces waves 
and geometric stability; we find some explicit results when we specialize the 
analysis to a Mooney-Rivlin solid with Newtonian viscosity. 
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2 Basic equations 



2.1 Kinematics 

Consider a continuous body whose stress-free reference configuration is denoted 
B r and in which material points are labelled in terms of their position vectors X. 
In the current (deformed) configuration at time t, denoted B, X occupies the 
position x, and the motion from B r to B is described by the bijection mapping 
X such that 

x = X (X,t). (1) 
The deformation gradient associated with the motion, denoted F, is defined 

as 

F = Gradz, (2) 

where Grad is the gradient operator in B r , and the velocity v of a material 
particle is defined as 

It follows that dF/dt = LF, where 

L = gradu, (4) 

with v regarded as a function of x and t, is the velocity gradient. Its symmetric 
part is the strain-rate tensor D, given by 

D = \{L + L T ), (5) 

where the superscript T denotes the transpose. Finally, the left and right 
Cauchy-Green deformation tensors are defined by 

B = FF T , C = F T F, (6) 

respectively, and we note that 

dC/dt = 2F T DF. (7) 

For an incompressible material only isochoric motions are permitted, in 
which case the constraints 

det F = l, trL = trD = 0, (8) 

are enforced at all times. The latter condition is equivalent to 

divv = 0, (9) 

where div is the divergence operator in B. 

For an incompressible material we now define eight independent invariants 
of the two tensors B and D by 

h = ixB, h = tr^- 1 ), h = tr(DB), J 6 = tr (DB 2 ), 

7 7 = tr(D 2 ), 7 8 =tr(D 2 B), J 9 = tr (D 2 B 2 ), I w =tv(D 3 ), (10) 

noting that the invariants -Z3 = det B = 1 and I4 = trD = have been omitted 
from the list by virtue of ([6]) and (JSj) . 
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2.2 Constitutive law and equation of motion 

For an incompressible isotropic material with Cauchy stress tensor T depending 
on B and D only, the general representation for the constitutive law is [T51 
Chap. 11] 

T = -pi + aiB + a 2 B 2 + a 3 D 

+ a 4 (DB + BD) + a 5 {DB 2 + B 2 D) + a 6 D 2 

+ a 7 (D 2 B + BD 2 )+a s (D 2 B 2 + B 2 D 2 ), (11) 

where J is the identity tensor, p is the Lagrange multiplier associated with the 
incompressibility constraint, and aj, i € {1, 2, . . . 8}, are material functions that 
depend on the eight invariants (|10[) : 

cti = a>i(Ii, h, h, h, I7, h, h, ho)- (12) 

The equation of motion in the absence of body forces is 

divT = pd 2 x/dt 2 , (13) 

where p is the mass density of the material. Equivalently, it can be written as 

DivS = pd 2 x/dt 2 , (14) 

where Div is the divergence operator in B r , and S is the nominal stress tensor, 
defined here as 

S = F- 1 T. (15) 



2.3 Equilibrium 

Suppose now that the material is in equilibrium in a deformed configuration B 
so that v = 0, D = 0. Let all quantities associated with B be denoted by an 
overbar. Then the deformation is written 

x = x(X), (16) 

the associated deformation gradient is F, and the corresponding left Cauchy- 
Green tensor is denoted B. The Cauchy stress is 

T= -pi + a!B + a 2 B 2 , (17) 

where 

a t =at i (h,l2,0,...,0), (18) 

and /1, J2 are the first two principal invariants of B. Finally, the equilibrium 
equation may be written in either of the equivalent forms 

divT = 0, DivS = 0, (19) 

where S = F~ 1 T. 
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3 Small motion superimposed on a static finite 
strain 



3.1 Incremental kinematics 

We now superimpose a small amplitude motion on the finite static deformation 
in the configuration B. Let x(X,t) denote this incremental motion. We then 
change variables from (X, t) to (x, t) and introduce the mechanical displacement 
vector u(x,t) defined by 

x(X,t)=u(x(X),t), (20) 

there being no need to distinguish between x and x in this linearization. The 
corresponding increment in the deformation gradient is 

F = Grada; = HF, (21) 

where H = gradtt is the displacement gradient. 

Because the basic deformation is static, we have v = d(x + x)/dt = du/dt, 
and hence 

dF/dt = Gr&dv = LF, (22) 
where L = gradu is the velocity gradient defined in (U), and we have 

L = dH/dt. (23) 

We compute the (linearized) increments in the relevant kinematical quanti- 
ties as 

B = HB + BH T , (B 1 ) = HB 2 + BHB + BH T B + B 2 H T , 
ix=2tx{HB), i 2 = -2ti (HB^ 1 ), 

h = tr (DB), i 6 = tr (DB 2 ). (24) 

Note that the increments in the invariants It, . . . , iio are zero at first order. In 
fact, because D = (L + L T )/2 is infinitesimal by (|23p . we have 

I 7 = I s =I g = I w = (25) 

at first order. It follows that the material parameters ai in the constitutive 
equation (jlip need from now on be considered as functions of four invariants 
only, namely (Ji, 7 2 , h, Thus, 

a i = a i (I 1 ,I 2 , h, h), ie {!,..., 8}. (26) 
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3.2 Incremental stress and incremental equations of mo- 
tion 

Now we increment the constitutive law (|lip . retaining only the first-order terms. 
We find, using (|24[) and the increment of (|26p . that 

T = - pi + fii (HB + BH T ) + a 2 (HB 2 + BHB + BH T B + B 2 H T ) 

+ a 3 D + a 4 (DB + BD) + a 5 (DB 2 + B 2 D) 

+ [2a n tr (HB) - 2a X2 tr (H B^ 1 ) + a 15 tr{DB) + a 16 ty(DB 2 )]B 

+ [2a3itr(HB) - 2a 22 tr(i? J B _1 ) + a 25 tr{DB) + a 26 tr(DB 2 )]B 2 , 

(27) 

where on and atj are the values of ai and don /dlj , respectively, evaluated for 
B = B, D = 0. 

Incrementing the connection (115[) between Cauchy stress and nominal stress, 
and using (|21|) , we obtain the increment in the nominal stress as 

S = F~ 1 (f-Hf). (28) 
It follows that the increment in the equation of motion (| 14[) , which is 

DivS = P d 2 x/dt 2 , (29) 

can equivalently be written as 

div (f - HT) = pd 2 u/dt 2 , (30) 

with x and t as the independent variables. This is coupled with the incremental 
incomprcssibility condition 

divit = 0. (31) 

We recall that the underlying deformation is homogeneous so that B is 
constant, and hence on and a^- are constants, as is T . It follows that 

div(JIT) = r(divH) = Tgrad(divtt) = 0, (32) 

and the equation of motion ([50)) reduces to 

div f = pd 2 u/dt 2 . (33) 

Upon substitution of (I27p into (f3"3"l) , we arrive at 

- gradp + div [aiBH T + a 2 (BHB + BH T B + B 2 H T )j 

+ div[a 3 D + a 4 {DB + BD) + a 5 (DB 2 + B 2 D)} 

+ Bgrad[2antr(i?B) - 2a 12 tr(HB~ 1 ) +a 15 tr(DB) + a 16 tr(DB 2 )} 

+ B 2 grad [2a 2 itr (HB) - 2a 22 tr (HB 1 ) + a 25 tr (DB) + a 26 tr (DB 2 )] 

= pd 2 u/dt 2 , (34) 

where we have used the incremental incompressibility condition IpTl]) . 
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3.3 Two-dimensional waves 

Let Af , A|, A| be the eigenvalues of B and denote by Xi,x 2 ,x 3 the coordinates 
associated with Cartesian axes along the corresponding eigenvectors. These are 
the principal axes of pre-strain, and for isotropic materials, as considered here, 
they are aligned with the principal axes of the pre-stress. 

In the remainder of the paper we focus on two-dimensional waves, whose 
spatial variations depend on two principal space variables only, x\ and X2 say. 
Hence 

u = u(xi,X2,t), p = p(xi,X2,t), (35) 
and the incremental incompressibility constraints (|3 1 1) and © reduce to 

"1,1 + "2,2 = 0, +"2,2 = 0, (36) 

respectively, where the comma denotes partial differentiation. The components 
of T in the {x\, X2) plane are 

f n =-p + 2(ai + 2A?a 2 )A?ui,i + (a 3 + 2\\a A + 2Aia 5 )"i,i 
+ 2A?(a u + A?a 2 i)(A?ui,i + A^ 2 , 2 ) 

- 2Af (ai2 + A?a22)(A7/ 2 tti,i + X^u^^) 
+ A?(ai5 + A?q:25)(A?wi,i + \\v 2 ,2) 

+ A?(ai6 + Aia 2 6)(Af«i,i + ^2,2), 
T22 = - p + 2(ai + 2A2a 2 )A2U 2 ,2 + (S3 + 2A^a 4 + 2A2<5 5 )w2,2 
+ 2A|(an + Ala 2 i)(Aftti,i + A|u 2 ,2) 

- 2Al(ai2 + A2a 22 )(A ] " 2 ui4 + X 2 2 u 2 ,2) 
+ A|(ai 5 + A|a 2 5)(Aii;i,i + A^2,2) 

+ >2(ai6 + A|a 2 6)(Afwi,i + A|u 2 , 2 ), 
T12 = [on + (A? + A|)a 2 ]A|ui, 2 + [ai + (A? + A^)a 2 ]w 2 ,i 

+ I [as + (A? + Ai)a 4 + (Af + Ai)a 5 ]("i,2 + "2,1), (37) 

and they do not involve U3. 

The incremental equations of motion ()33|) reduce to 

Tii,i + Ti2, 2 = P"i,tt, Ti 2l i + T 22 ,2 = pu 2 ,tt- (38) 

It is easy to check that these equations decouple from the third equation of 
motion Ti3 : i+T 2 3,2 = P"3,tti which involves U3 only. We therefore take 7/3 = so 
that this is satisfied and need not be considered further. A simple manipulation 
of (f3"8")) then leads to 

(Tn - T22X12 + T12, 22 - Ti2,n = p(ui,2tt - U 2 ,ltt), (39) 
which eliminates p. 
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It is now convenient to introduce the material parameters a, 7, /3, 6, e, defined 

by 

a = [&i + a 2 {\\ +X\)\\\, 
1 =[a 1 +a 2 {\\ + \l)]\l, 
2/3 = ai{\\ + A|) + a 2 (3\f + 3A^ - 2A 2 A 2 ) 

+ 2<5n(A 2 — A 2 ) 2 + 2ai 2 A^ 2 A2~ 2 (A 2 — A 2 ) 2 
+ 2a 21 (A 2 - A 2 ) 2 (A 2 + A 2 ) + 25 22 Ar 2 A 2 - 2 (A 2 - A 2 ) 2 (A 2 + A 2 ), 
25 = a 3 + ai{\\ + A 2 ) + a 5 (Xf + A 2 ), 
e = [a 15 + (a 16 + a 25 )(A 2 + A 2 ) + a 26 (A 2 + A 2 ) 2 ](A 2 - A 2 ) 2 . (40) 

Then, on use of u 2 , 2 = -fii.i and v 2 _ 2 = — «i,i, we obtain 

Tn - T 22 = {a + 7 + 20)ui tl + (e + 45)^,!, (41) 
T12 = aii2,i + 71*1,2 + 5(«i,2 + V2,i). (42) 

The incremental incompressibility constraint suggests the introduction of a 
scalar potential function ip = ip(xx, x 2 , t) such that 

wi =- 0,2, u 2 = -ip,i, Vi=ip t 2t, v 2 = -ip,iu (43) 

use of which enables the equation of motion (|39p to be cast as an equation for 
if), namely 

Chilli + 2 ) 9V,1122 + 7^,2222 

+ <KV>,llHt + 2^,1122* + V\2222t) + £"0,1122* = P(lp,Utt + lp ,22tt) ■ (44) 

This is the equation that governs the two-dimensional incremental motions. 
3.4 Dissipation 

For a continuum, the work done by external forces is converted into kinetic 
energy, stored energy, and dissipated energy. The combination of the latter two 
is measured by the rate of working of the stresses, which, per unit volume, is 
tt(SF). For the incremental motion this can be written tr[(S + S)F]. The 
first term in this sum can be considered as the stored elastic energy associated 
with the underlying static deformation, whilst the second term tr (SF) is a 
measure of the dissipation associated with the motion (which may include some 
additional stored energy). 

From (|28p . ([22]) . the symmetry of T, and the definition ([5]), we obtain 

tr(SF) = tr(TD) -tr (HTL). (45) 

For the two-dimensional incremental motions, the two terms on the right-hand 
side of (|45|) may be computed, respectively, as 

tr (TD) = T n D u + T 22 D 22 + 2T 12 D 12 

= (Tn ~ T 22 )vi t i + Ti 2 (ui,2 + V2,i) 
= (a + 7 + 2/3)ui ) i« u + (e + 4c5)u 2 a 

+ (cra 2: i + 7^1,2X^1,2 + v 2 ,i) + <5(vi, 2 + w 2 ,i) 2 (46) 
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and 

tx(HTL) = Th(mi,iWi,i + u 2 , 1^1,2) + ^(^i, 2^2,1 + u 2 . 2 v 2 , 2 ) 

+ ri 2 (ui,2Wl,l + 1*2,2^1,2 + Ul,l«2,l + U2,l«2,2)- (47) 

In the case of time-periodic homogeneous motions, we use angle brackets to 
denote the time average over a period; here we find that 

(tr(TX>)) = (e + 45)(vl 1 )+S((v h2 + v 2A ) 2 ), (48) 

and that the other terms vanish, as we now show. 
First we have 

Ul, = ip,i2ip,i2t = \ (V 2 i 2 ) 1 > ( 49 ) 
whose time average clearly vanishes by periodicity. Next we have 

(cm 2 ,i+ 7^1,2) («l,2+«2,l) = 5«('/' 2 li) t + \l (3^22) t -aV',iiV',22t-7 ? /',22V',iit- 

(50) 

Here the time averages of the first two terms vanish by periodicity. To com- 
pute the time averages of the last two terms, we write ip explicitly. For (two- 
dimensional) time-harmonic motions, we may write it in the general form 

ip = Cc iuj{s n - x - t] + rje- 1 "^"-" 5 -*) , (51) 

where C is a complex constant, uj is the real frequency, s is the complex slowness, 
n = (m, n 2 , 0) is a real unit vector in the propagation direction, and the overbar 
denotes the complex conjugate. Introducing the function ip defined by 

(52) 

we obtain the expressions 

^,11^,22* = ip,n^,22t = \n\n\ (ip 2 ) t , (53) 

which have a zero time average by periodicity. Similar calculations show that 
the time average of tr (HTL) vanishes. 

Turning back to the time average of (|4"8"|) , we find, using the function ip, that 
it can be written as 

(tv(fD)) = (en 2 1 nl + 5)(( lp ^) 2 ), (54) 

making it clear that the deformed viscoelastic solid is dissipative under (plane) 
incremental motions when en 2 n 2 + S > for all m, n 2 such that n\ + n 2 , = 1. 
This is ensured when 

6>0, e + 4<S>0, (55) 

with at least one of these inequalities being strict. For the remainder of the 
paper, we assume that these inequalities hold. 

Having established the conditions for time-averaged dissipation of time- 
periodic homogeneous motions, we now investigate stability issues for a de- 
formed viscoelastic solid occupying, first, the whole space and, second, a semi- 
infinite space. 
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4 Material stability 



First we look at the situation where the perturbation has no time dependence, 
that is when d/dt = 0. For all intents and purposes, the viscous effects are not 
felt then, and the solid behaves as a purely elastic solid. The corresponding 
incremental equation of elastostatics is the specialization of to 

a^ini + 2/3^1122 + 7^,2222 = 0. (56) 

It is known [11] that this equation is strongly elliptic when 

a > 0, 7 > 0, (3 + y/dn>0, (57) 

and we assume henceforth that these inequalities hold. This guarantees ma- 
terial stability in the strong ellipticity sense with respect to incremental static 
deformations. 

Next we study bulk homogeneous plane waves because they provide a natural 
tool for addressing the question of the material (bulk) stability of the deformed 
viscoelastic solid. We therefore seek solutions of the form 

iP^iPoe^ 3 ™-^, (58) 

where tJjq is a constant, lu = lu + + ilu~ is the complex frequency, s = s + + is - 
is the complex scalar slowness, and n is a real two-dimensional unit vector in 
the direction of propagation. Note that this motion is not necessarily time- 
periodic because lu~ may be different from zero. Combining this motion with 
the expressions in (|43|) . we see that the displacement, velocity and stress fields 
have the same exponential dependence as ip. The argument of the exponential 
may be decomposed as 

iu)(s n • x — t) = — [(cj + s~ + uj~s + ) n • x — uu~t\ 

+ i [(w + s + — oj~s~) n • x — lo + {\ . 

The first bracketed term gives the amplitude variations of the fields, and the 
second one their phase. 

Material stability is ensured when there is no amplitude growth for a given 
phase [5]. In other words, when (uj + s~ + u~ s + ) n • x — uj~t > with (lu + s + — 
oj~s~)n • x — u + t = 0. This gives 

"/ S ; W V -c->0, (59) 
uj + s + — lu s 

or cquivalently, after removing the positive factor [(^ ,+ ) 2 + (u; - ) 2 ], and taking 
the inverse, 

s+ 

—lu + -lj->0. (60) 
s 

We now examine the implications of this inequality for a deformed viscoelastic 
solid. 
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Substitute the expression (|5"5j) for ijj into the equation of motion ([55]) . and 
separate the real and imaginary parts to obtain 

anf + 2(3n\n 2 2 + -ynj + uj~(5 + en\n\) = p[(v+) 2 - (u~) 2 ], 

uj+(5 + enlnl) = -2pv+ir, (61) 

where are real quantities defined by v + + iv~ — (s + + is - ) -1 , that is 

(62) 



(s+) 2 + (s-) 2 ' (s+) 2 + (s-) 2 ' 

From equation (|6ip ? we find 

v~ u> + (6 + en\n 2 ) 
~ ~~ 2 /9 (w+) 2 



(63) 



Then, dividing equation (f6T]) i through by (u + ) 2 , and using this latter identity 
we find an expression for lj~ . We can also use the identity above to find s + / s~ = 
—v + /v~. We end up with 



p(v + ) 2 — anf — 2f3n\n\ — jn 2 1 
5 + en 2 n 2 



--oiv+fiu+fiS + enlni) 



£_ = 2 _Piv_)_ 

s~~ 5 + en\n 2 



2 

2 



^ + ~ = 2 x ?,X 2 - ( 64 ) 



Hence the stability condition (|6H)) reads 

p(v + ) 2 + anf + 23n 2 n 2 , + 777,9 , ,. 0/ 1.0,,. 9 9s / 

2 r 2 + p(t; + ) 2 (cj + ) 2 (<y + enfn 2 ) > 0. (65) 



This condition is clearly satisfied when both (|55[) and (|57|) hold. In other 
words, time-averaged dissipation with respect to time-periodic motions, coupled 
to strong ellipticity with respect to static deformations, results in material sta- 
bility. 

Before we go on to investigate geometric stability, we pause to consider a 
classic sub-case of the general bulk wave (158[) . namely the damped travelling 
wave solution. It is of the form 

ip = </;oe~ at cos k(n ■ x - ct), (66) 

where a is the damping factor, k is the wavenumber, and c is the speed. It is 
an important subclass of (|58|) . obtained by taking 

u + s- +uj- S + = 0, (67) 



so that there is no spatial attenuation of the amplitude. Then |58J specializes 
to (|66|) by making the identifications 

a = — uj~ , k = —a/v~, c = v + . (68) 
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Also, (|67[) gives cu + = (s + / s )a = —(v + /v )a, so that the dispersion equations 
(|6"Tj) reduce to 

pc 2 = an\ + 2fin\n\ + -(5 + en\n 2 2 ) 2 'k 2 / '(4p) , 

2pa = (5 + enjn^k 2 . (69) 

If damped travelling waves (f66|) can be generated in a viscoelastic solid, 
then these equations provide a means to determine the constitutive parameters 
by variation of the propagation direction and of the underlying deformation. 
In particular, if the response of the solid shows a dependence of the damping 
factor a on the direction of propagation, then the constitutive model must be 
such that e ^ 0, according to (|69|) 9. By (|40)) .^. this means that the constitutive 
parameters a.\ and cannot both be completely independent of the invariants 
is and Iq. Therefore, only certain (quite complex) constitutive models can 
display an influence of the propagation direction on the damping factor. If the 
model is such that e > 0, then the directions of maximal dissipation are along the 
bisectors of the principal directions, and those of minimal damping are aligned 
with the principal axes (and vice- versa if the model is such that e < 0). 

Conversely, if the response of the solid shows that the damping factor a is 
independent of the direction of propagation for damped travelling waves, then 
the constitutive parameters ot\ and can both be completely independent of 
the invariants I5 and I§. 

5 Geometric stability 

To study surface stability, we consider inhomogeneous motions in the half-space 
X2 > with boundary x 2 — in the (x%, a^-plane, the deformation correspond- 
ing to pure homogeneous strain with the principal axes of strain coincident 
with the Cartesian axes. On the surface X2 = we assume that the incremental 
surface tractions vanish, i.e. 

{T - HT) 21 = (T — HT) 22 = 0. (70) 

The shear traction condition leads, after some manipulation, to a condition 
involving -0, namely 

7(^,22 - + 0-2*0,11 + <W,22i - 1p,llt) = (71) 

on 12 = 0, where a 2 is the (uniform) principal stress normal to the boundary 
in basic state of deformation, i.e. <r 2 = T 22 . 

In order to express the normal component of the incremental traction on the 
boundary in terms of it is first necessary to differentiate the latter equation 
in (|70p along the boundary, i.e. with respect to xi, and then make use of the 
first component of the equation of motion to eliminate dp/dxi (assuming this 
holds on the boundary). After further manipulations this leads to 

(2/3 + 7 - cr 2 )V\ii2 + 7^,222 + (e + 3c5)^,ii2t + -5-0,222* - pip,2tt = (72) 
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on X2 = 0. When the viscous terms are absent (S = e = 0) the boundary 
conditions are then precisely those given by Dowaikh and Ogden [IT] for the 
purely elastic case. 

We now consider waves of the inhomogeneous form 

t/j = ip e^ kxi - wt h- ksX2 , (73) 

where k, lu and s may be complex. However, we impose the following propagation 
inequalities 

Re(fc) > 0, Im(fc) > 0, Re(w) > 0, (74) 

so that the wave propagates in the positive x\ direction at the interface X2 — 
and attenuates (if at all) in the positive X\ direction. Additionally, we set 

Re(ks) > 0, (75) 

so that the wave decays away from the boundary x-i = (the localization con- 
dition). Finally, we pay special attention to the sign of Im(cj); clearly, if 

Im(cj) < 0, (76) 

then the wave is damped (decays in time); if Im(w) > it blows up in time, 
indicating the onset on instability, at least in the linearized theory. We refer to 
(|76p as the stability condition. If Im(w) = there is neither growth nor decay 
in time. 

On substitution of ([73")) into equation (|4"4")l we obtain a bi-quadratic for s, 
which can be written compactly in the form 

js 4 -(2(3-n)s 2 + a-n = 0, (77) 

where we have introduced the notations 

a = a-iwS, 2/3 = 2(3 - 2iw5 - iwe, 7 = 7-^, £l = puj 2 /k 2 . (78) 

The general solution of the equation of motion may be written in the form 

ip = e Wxi-"t)(Ae- ksiX2 + Be~ kS2X2 ), (79) 

where A and B are constants and s\ and S2 are the solutions of ([77)) that satisfy 
([74]) , (J75"j) . and Substitution of (JTSJ) into the boundary conditions fl7TJ) and 
(172^) then gives the two equations 

fi(sl + 1) - a 2 ]A + + 1) - <r 2 ]B = 0, 

[2/3 + 7 - a 2 - Cl - jsf]siA + [2/3 + 7 - a 2 - Cl - jsj^B = 0, (80) 
for A and B. 

After removal of the factor si — S2 the determinant of coefficients yields the 
dispersion equation, which, on use of the sum and product of the roots of (|77[) , 
reads 

(7 " ^2) 2 - 7(a -Cl)- 7Sis 2 (2/3 + 2 7 - 2a 2 - Cl) = 0. (81) 
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The product s\S 2 has not been replaced since there are two possible solutions 
of s\s\ = (d — f2)/7 and this needs careful evaluation. This dispersion equation 
reduces to the elasticity result obtained by Dowaikh and Ogden [TT] on setting 
(5 = e = and ui and k real. We remark here that if the case s\ = s 2 is consid- 
ered separately and the solution (|75|) amended accordingly then the associated 
dispersion condition reduces to a 2 = ±27. It can be shown that this also follows 
from the appropriate specialization of (|8Tj) . However, since a 2 is real and, in 
general, 7 is complex, this cannot be satisfied unless Re(o;) = or S = 0. In the 
purely elastic case the corresponding special solution is a 2 = ±27 [TT] . 
From (J75J) it follows that 

2/3-d-7 = 2/3-cv-7 - ieco. (82) 

Now, in the context of elasticity, materials for which 2/3 — a — 7 = form a 
special class and lead to simplifications in the analysis. Similar simplifications 
occur here if we focus on viscoelastic solids for which 



2/3 = a + 7, e = 0, (83) 

and we assume henceforth that the material model is specialized in accordance 
with (|83)). Then, {77]) factorizes to give 

(s 2 - l)[js 2 - (a - (l)} =0. (84) 

One root consistent with the restrictions ([74")l . (|75|) is s — 1 and we refer to this 
as s%. 

There are two possibilities for the second root, 




s 2 =±\— r - . (85) 



A test must be conducted by computing ks 2 for each possibility and checking 
whether the localization requirement (|75p is satisfied. Here the square root 
symbol designates the complex number with square equal to (d — fi)/7 and 
positive real part. In any case, the dispersion equation (|8ip can be re-cast as a 
cubic in s 2 , namely 

s^ + s^ + (3-2a 2 )s 2 -(l-<T 2 ) 2 = 0, (86) 

where a 2 = 02/7. Note that this cubic is not obtained by a squaring process, in 
contrast to the cubic obtained by Currie et al. [6] . It does not contain spurious 
roots a priori and it is therefore legitimate to check the validity of each of 
its three roots against conditions (f74|) . ([75]) . and (f76|) once k (or ui) has been 
deduced from (f85|) for a given uj (or k). 

However, the behaviour of the roots is highly dependent on the material 
parameters and on the pre-stress and pre-strain, and little can be concluded 
in general. In order to make progress and provide an illustrative example, we 
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first specialize the analysis further to a Mooney-Rivlin solid with Newtonian 
viscosity, with constitutive equation 



T = -pI+{C l +C 2 I 1 )B-C 2 B 2 + vD, (87) 

where C\, C%, and v are positive constants. Then the parameters a, 7, (3, S, e 
of (J3QJ) reduce to 

a = (Ci + C 2 X 2 3 )Xl, ~f = (d + C 2 X 2 3 )X 2 2 , 2/3= (Ci +C 2 A 2 )(A 2 +A 2 ), 
(5 = 0, e = 0, (88) 

making it clear that this solid belongs to the class (|83|) . The quantity fi = Ci+C 2 
is its static shear modulus and v is its dynamic viscosity. Next, we specialize 
the pre-deformation and pre-stress to a plane strain with no normal load, 

Ai = A, A 2 = A-\ A 3 = l, (T2-0, (89) 

where A is the stretch ratio in the x± direction. By taking v = 0, this set- 
up allows for direct comparison with the purely elastic Mooney-Rivlin case, for 
which Biot [3] showed that the critical compressive stretch for surface instability 
is A cr = 0.5437. Also, s 2 is now a root of the cubic s\ + s\ + 3s 2 — 1 = 0, 
independent of the material parameters and the pre-deformation. Explicitly, s 2 
is among the three solutions of this cubic, which are 

s° = 0.2956, s£ = -0.6478 ± 1.721 i, (90) 

and the dispersion equation is deduced from (|85| as 



/ A 2 — imv/ /i — puj 2 / (^k 2 ) 
S2 = ±\ r-2 : 7 ■ (91) 

Despite the strong simplifying assumptions made here, the possibilities for 
solutions to the surface stability problem remain rich and varied because of the 
possible complex nature of the wave number and of the frequency. 

5.1 Real frequency, complex wave number 

First we take u> real (co > 0). Then there is neither growth nor decay in time. 
In other words, taking u> real is not appropriate for the study of stability. Nev- 
ertheless, we may investigate the possibility of a surface wave existing, i.e. a 
solution in the form ([75)1 satisfying all four conditions (TT4")) . (|75p. When we 
choose s 2 = s§ = 0.2956 as the root from (|9T))) , we find that these conditions are 
equivalent to 

Bet [E!L) >0} ha(^-)>0. (92) 

VV puJ VV p u ) 

Figure [T] shows the variations of these quantities with respect to A for several 
values of the dimensionless parameter voj / fi. In the first figure, the dashed curve 



15 





vu/H = 




A 



Figure 1: Dimensionless slowness and damping factor for a surface wave with 
real frequency in a Mooney-Rivlin viscoelastic solid subject to plane strain. The 
analysis shows that the half-space becomes unstable when compressed by more 
than about 46% (vertical asymptote) and this confines the validity of the curves 
to the right of the vertical dashed line. 

represents the (non-dimensional) slowness in the purely elastic case {v = 0), with 
a vertical asymptote at the critical stretch A cr = 0.5436 where the speed drops 
to zero. The introduction of viscosity removes this singularity and a surface 
wave may propagate for the whole compressive range, unless of course the half- 
space becomes unstable (see below). In the second figure there is no curve at 
v = because the purely elastic wave is not damped [13] . 

When we take either si = as the root from (f9"0"]) . we find that the con- 
ditions (|74f . (|75|) cannot be satisfied simultaneously. Hence, there is only one 
possibility for a surface wave to propagate over a deformed viscoelastic Mooney- 
Rivlin solid, that which tends to the Rayleigh surface wave solution when the 
viscosity tends to zero. This is in accord with the results of Romeo [26] in linear 
elasticity (no finite pre-deformation) . 

5.2 Complex frequency, complex wave number 

When we allow both the frequency and the wave number to be complex, we find 
that the imaginary part of ui is negative only in the range where A 2 — A _2 s| > 
0, for S2 = s§ = 0.2956; the other two roots of ([90l do not yield any 
conclusion with respect to stability analysis. When A 2 — A~ 2 s 2 < 0, i.e. when 
A < A cr = 0.5436, the imaginary part of u> is positive, indicating instability. 
When A = A cr — 0.5436, both the real and imaginary parts of uj are zero, as can 
be checked from the dispersion equation (|9"T]) . The conclusion is then that the 
half-space becomes unstable when the complex speed uj/k is zero, just as in the 
purely elastic case. This is in accord with the correspondence principle of Biot 
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